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A comparative study is made on the metal-insulator transition of Dirac fermions in the honeycomb and tt- 
flux Hubbard models at half filling by means of the variational cluster approximation and cluster dynamical 
impurity approximation. Paying particular attention to the choice of the geometry of solver clusters and 
the inclusion of particle-bath sites, we show that the direct transition from the Dirac semimetallic state to 
the antiferromagnetic Mott insulator state occurs in these models, and therefore, the spin liquid phase is 
absent in the intermediate region, in agreement with recent quantum-Monte-Carlo-based calculations. 


1. Introduction 

The mechanism of the metal-insulator transition 
(MIT) has long been one of the central issues in strongly 
correlated electron systems.^’In particular, the MIT 
in correlated Dirac fermion systems has attracted much 
attention recently, a typical example of which is the 
honeycomb-lattice Hubbard model at half filling repre¬ 
senting graphene.Because the honeycomb lattice is bi¬ 
partite and free from frustration, the Neel antiferromag¬ 
netic (AF) Mott insulator (MI) state is realized in the 
strong coupling region. However, unlike in the square- 
lattice Hubbard model, where perfect Fermi surface nest¬ 
ing is present, one expects that the AF order will not ap¬ 
pear in the weak coupling region but rather the massless 
Dirac semimetallic (SM) state will be maintained until a 
critical interaction strength is reached. 

The MIT in the honeycomb lattice was studied by 
Meng et al.^^ using the quantum Monte Carlo (QMC) 
method, whereby they claimed the presence of a quan¬ 
tum spin liquid (SL) state (or nonmagnetic MI state) 
in the intermediate region between the Dirac SM state 
and the antiferromagnetic Mott insulator (AFMI) state. 
Their study attracted much interest because it suggested 
the emergence of the SL state in systems without frus¬ 
tration in their spin degrees of freedom. However, subse¬ 
quent studies based on the large-scale QMC method, 
the pinning field approach using the QMC method, 
and analysis of the quantum criticality by finite-size scal- 
ing^^’ have consistently suggested the direct transition 
from the SM state to the AFMI state, and therefore we 
now anticipate that the SL state is absent in this model. 
Similar debates have also been had for the 7r-flux Hub¬ 
bard model, another Dirac fermion system, whereby the 
direct transition from the SM state to the AFMI state is 
now anticipated.^^ 

Quantum cluster methods have also been used to study 
the MIT in the honeycomb Hubbard model.In par¬ 
ticular, cluster dynamical mean-field theory (CDMFT) 


and variational cluster approximation (VCA) calcula¬ 
tions have shown that if the 6-site hexagonal ring is 
used as a solver cluster, the single-particle band gap 
opens even in the weak coupling region where the AF or¬ 
der is absent, thereby suggesting the presence of the SL 
state.However, the opening of the band gap at the 
infinitesimal interaction strength was questioned, 
and moreover, from comparison with the results of the 
cluster dynamical impurity approximation (CDIA) and 
dynamical cluster approximation (DCA), the emergence 
of the nonmagnetic insulator phase predicted by the 
CDMFT and VCA was considered to be unrealist 
So far, not much is known about the MIT in the 7r-flux 
Hubbard model studied by quantum cluster methods. 

In this paper, motivated by the above development in 
the field, we will make a comparative study on the MIT 
of correlated Dirac fermions in the honeycomb and 7r-flux 
Hubbard models at half filling by means of the VCA and 
CDIA. We will, in particular, point out that a suitable 
choice of the cluster geometry is essential in the quantum 
cluster calculations to suppress the opening of the band 
gap in the weak-coupling region and that the inclusion of 
particle-bath sites is important in discussing the order of 
the MIT as well as the transfer of spectral weight in the 
single-particle spectral function. We will thereby show 
that the direct transition from the Dirac SM state to the 
AFMI state occurs with increasing interaction strength 
in these models and that the SL phase is absent in their 
intermediate coupling region. 

2. Models and Methods 

The honeycomb and tt-Hux Hubbard models may be 
defined by the Hamiltonian 

- X T Cja + t/ X > (1) 

i 

where c\^ is the creation operator of a fermion (which will 
be referred to as an electron hereafter) with spin a at site 
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Fig. 1. (Color online) Schematic representations of the (a) hon¬ 
eycomb it' jt = 0) and (b) 7r-flux it' jt = —1) lattices. The dashed 
line in (a) indicates the bonds with the hopping parameter t' and 
the red lines in (b) indicate the bonds with the negative hopping 
parameter t' = —t. Noninteracting DOSs [(c) and (d)] and contour 
plots of the band dispersions [(e) and (f)] are also shown for 
the honeycomb (left panels) and 7r-flux (right panels) lattices. The 
green dots in (e) and (f) indicate the Dirac points in fc-space. 


Fig. 2. (Color online) Single-particle gap Asp in the PM state 
and local magnetization m in the AF state calculated by the VC A 
as functions oiU/t. The results for the honeycomb Hubbard model 
are shown in (a), (b), and (c) and those for the 7r-flux Hubbard 
model are shown in (d), (e), and (f). The geometry of the solver 
cluster used is illustrated in each panel. 


i and nia = c\^Cicj. Uj is the hopping amplitude: we define 
tij = t for the nearest-neighbor bonds and Uj = t' for the 
bonds connecting hexagons in the honeycomb lattice [see 
Fig. 1(a)]. U is the on-site Coulomb repulsion. We assume 
the filling of one electron per site (half filling). Changing 
the value of t' , we can tune the system continuously from 
the honeycomb lattice at t' = 0 to the 7r-fiux lattice 
at t' = —t [see Fig. l(b)].^^^ At U = 0, these systems 
at low energies are described in terms of the massless 
Dirac fermions; their densities of states (DOSs) and band 
dispersions are shown in Figs. l(c)-l(f). 

We apply the which is a quantum clus¬ 

ter method based on self-energy functional theory 

(SFT).29^30) 

we introduce disconnected 
finite-size clusters (that are solved exactly) as a refer¬ 
ence system. By restricting the trial self-energy to the 
self-energy of the reference system S', we can obtain the 
grand potential of the original system in the thermody¬ 
namic limit as 

n = n' ^Tt ln(Go ^ - S')-^ - Tr In(G'), (2) 

where G', and Go are the grand potential, the Green’s 
function of the reference system, and the noninteract¬ 


ing Green’s function, respectively. The short-range cor¬ 
relations within the cluster of the reference system are 
taken into account exactly. The one-body parameters t' 
of the reference system are optimized according to the 
variational principle {t^)]/dt' = 0. In the VGA, we 

can treat the spontaneous symmetry breaking by adding 
appropriate Weiss fields to the reference system.We 
have to choose an exactly solvable reference system; here 
we apply an exact diagonalization method and solve the 
quantum many-body problem in the cluster of the refer¬ 
ence system. 

We also use the CDIA,^^^ which is an extended ver¬ 
sion of the VGA where particle-bath sites are added to 
the clusters to take into account the electron-number 
fiuctuations in the correlation sites. In the GDIA, we 
optimize the hybridization parameter between the bath 
and correlation sites V and the on-site energy of the 
bath sites 5 based on SFT.^^’^^^ Note that the GDIA 
is intrinsically equivalent to GDMFT with an exact- 
diagonalization solver. 

3. Results and Discussion 

3.1 Results of VC A 

First, let us consider the single-particle gap Agp in 
the paramagnetic (PM) state and the staggered mag- 
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Fig. 3. (Color online) Single-particle spectral function A(k^u) 
[(a) and (b)] and DOS N(u;) [(c) and (d)] in the PM states of 
the TT-flux Hubbard model calculated by the VGA with the 12-site 
cluster. The horizontal dashed line in (a) and (b) indicates the 
Fermi level. We applied the artificial Lorentzian broadening of the 
spectra of r]/t = 0.15 in (a) and (b) and r]/t = 0.05 in (c) and (d). 
The inset in the lower panels is an enlargement of the DOS near 
the Fermi level, assuming rjjt = 0.005. 


netization m in the AF state, which are calculated for 
the honeycomb and tt-Hux Hubbard models using the 
VGA. We introduce the Weiss field associated with the 
two-sublattice Neel order and evaluate the local magne¬ 
tization m = — Tiii). The gap Agp is evaluated in the 

absence of the Weiss field as the jump of the chemical 
potential with respect to the number of electrons in the 
system. Note that the band gap always opens when the 
AF order appears. We use clusters of 6, 10, and 12 sites 
as reference systems; the clusters used for the honeycomb 
and TT-flux lattices are topologically equivalent but with 
different hopping parameters (see Fig. 2). 

The results for the honeycomb Hubbard model are 
shown in Figs. 2(a)-2(c). We find that the MIT is sensi¬ 
tive to the choice of the clusters, i.e., the results obtained 
using the clusters of 6 and 12 sites are qualitatively dif¬ 
ferent from those in the case of 10 sites. The AF order 
appears at I/af/^ = 3.8 for the clusters of 6 and 12 sites, 
the results of which are in good agreement with results 
of QMC simulations.^’ However, the gap Agp opens at 
infinitesimal U values and the SM phase appears only at 
1/ = 0, thus suggesting the presence of the PM insulator 
state at 0 < 1/ < Uaf- Recent studies, however, have 
claimed that this gap cannot be regarded as the true 
Mott gap,^^’^^^ the details of which will be discussed be¬ 


low. For the cluster of 10 sites, on the other hand, the 
SM phase persists up to a large U value and the transi¬ 
tion to the AF phase occurs directly from the SM phase. 
Here, the AF order appears at 1 /af/^ = 2.7 and the gap 
Agp opens at Upu/t = 3.0, qualitatively consistent with 
the results of recent QMC simulations, where the direct 
transition from the Dirac SM phase to the AFMI phase 
was predicted.^’ 

The results for the TT-flux Hubbard model are shown in 
Figs. 2(d)-2(f). We find that the results obtained using 
the clusters of 6, 10, and 12 sites are qualitatively the 
same as each other, i.e., the SM phase persists up to a 
large U value. The AF order appears at I/af/^ = 3.4, 

2.9, and 3.4 and the gap Agp opens at Upu/t = 4.5, 

4.9, and 4.8 for the clusters of 6, 10, and 12 sites, re¬ 
spectively. Therefore, the PM insulator state does not 
exist between the SM and AFMI phases, in accordance 
with the results of recent QMC simulations that show 
the direct transition from the SM phase to the AFMI 
phase.Note that the transition point Uaf of our 
VGA calculations is smaller than that of the QMC sim¬ 
ulations, UAF/t = 5.25-5.5,^^^ which may be due to the 
anisotropy of the clusters used in our calculations; the 
agreement becomes good if we use an isotropic cluster of 
4 sites, which gives the value UAF/t = 5.0. 

We also calculate the single-particle spectral function 
and DOS using cluster perturbation theory (CPT)^^^ for 
the PM state of the TT-flux Hubbard model. The results 
are shown in Fig. 3. We immediately find that the Dirac 
linear band dispersion is clearly visible near the Fermi 
level at 1/ < Ufm [see Fig. 3(a)], whereas the band gap 
opens at 1/ > Ufm [see Fig. 3(b)]. The transfer of spectral 
weight occurring with increasing U/t is seen in Figs. 3(c) 
and 3(d), which is characteristic of the VGA and will 
be discussed below in Sect. 3.2 in comparison with the 
results of the CDIA. 

3.2 Results of CDIA 

Next, let us discuss the roles of bath sites in MIT using 
the CDIA. Following a previous study on the honeycomb 
Hubbard model,we examine the honeycomb and tt- 
flux Hubbard models using the 4-site 6-bath cluster. The 
results are shown in Fig. 4. We find that the grand poten¬ 
tials of the honeycomb and TT-flux Hubbard models both 
have two stationary points around the transition point 
I/c- The SM solution exists at a small I/, which vanishes 
at Uc 2 with increasing U. The MI solution exists at a 
large I/, which vanishes at Ud with decreasing U. The 
two solutions thus coexist in the region Ud ^ U < Uc 2 , 
and the ground-state energies cross at Uc. We obtain 
the values Ud/t = 6.6, Uc 2 /t = 7.7, and Uc/t = 7.5 
for the honeycomb Hubbard model and Ud/t = 8.6, 
Uc 2 /t = 10.3, and Uc/t = 9.8 for the tt-Hux Hubbard 
model. The calculated result for Agp [see Figs. 4(c) and 
4(d)] shows hysteresis between the SM (Agp = 0) and MI 
(Agp > 0) solutions, which indicates that Agp jumps dis- 
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Fig. 4. (Color online) Ground-state energies [(a) and (b)] and 
single-particle gaps [(c) and (d)] in the PM state of the honeycomb 
and TT-flux Hubbard models calculated by the CDIA, where the 
4-site 6-bath cluster shown in the left panels is used. 


continuously at Uc- These results clearly indicate that the 
MIT is of the first-order (or discontinuous) in the CDIA, 
which is in contrast to the results of the VGA where 
the second-order (or continuous) transition is found (see 
Fig. 2). The first-order MIT is thus expected in the ac¬ 
tual honeycomb and 7r-flux Hubbard models in which the 
electron-number fluctuation is present. 

To further clarify the roles of bath sites in the MIT, we 
examine the U dependence of the single-particle spectral 
function and the DOS calculated using CPT. The results 
for the TT-flux Hubbard model obtained in the CDIA are 
shown in Fig. 5, where the Dirac linear band dispersion is 
clearly visible in the vicinity of the Fermi level. Note that 
the slope of the dispersion at the Dirac point becomes 
steeper for larger values of U. 

Comparing the DOS curves, we find that the results 
in the VGA (see Fig. 3) are indeed significantly differ¬ 
ent from those in the CDIA (see Fig. 5) in the following 
respects, (i) The spectral weight with a large peak at 
cj/t = 2 in the U/t ^ 0 limit [see Fig. 1(d)] is partially 
transferred to a broad higher-energy region correspond¬ 
ing to the “upper Hubbard band” with increasing U/t, 
which is observed in both the VGA and CDIA. (ii) With 
increasing I/, the remaining spectral weight at uj/t 2 
shifts to higher energies in the VGA [see Figs. 3(c) and 
3(d)], while in the CDIA, it shifts rapidly to lower en¬ 
ergies and simultaneously loses its weight [see Figs. 5(c) 
and 5(d)]. (iii) We thus have a large spectral weight at 
low energies {uu/t < 1) in the CDIA, which is rather 
small in the VGA. The spectral weight characteristic of 
the massless Dirac SM dispersions can, however, be seen 
in the vicinity of the Fermi level in both the VGA and 



Fig. 5. (Color online) As in Fig. 3 but for the results of the 
CDIA with the 4-site 6-bath cluster. In (a)-(c), the single-particle 
spectral functions and DOSs of both the SM and MI states are 
shown at U/t = 9, while in (d), the DOS of only the SM state is 
shown. 


CDIA spectra, (iv) More quantitatively, a kink appears 
in the lowest-energy region of the DOS in the VGA (see 
the inset of the lower panels of Fig. 3), which shifts to¬ 
ward the Fermi level with increasing U. The DOS curve 
becomes steeper near the Fermi level (or the /^-linear dis¬ 
persion becomes flatter at the Dirac point), renormaliz¬ 
ing the Fermi velocity but keeping the electrons massless. 
No quasiparticle peak appears. At a critical U value, the 
kink disappears and simultaneously the gap begins to 
open gradually. In the CDIA, similar but stronger effects 
can be seen with increasing U/t \it the lowest-energy re¬ 
gion of the DOS, until the gap opens discontinuously at 
Uc [see Figs. 5(c) and 5(d)]. These low-energy behaviors 
in the CDIA are consistent with the results of the single¬ 
site DMFT for the honeycomb Hubbard model^T) ^j^d 
are expected to be realistic in the honeycomb and 7r-flux 
Hubbard models where the electron-number fluctuates. 

3.3 Cluster geometry dependenee 

Finally, let us discuss the cluster geometry depen¬ 
dence of the single-particle gap Agp in the PM phase. 
In Figs. 6(a) and 6(b), we show the results of the 6-site 
6-bath system for the honeycomb Hubbard model and of 
the 4-site 4-bath and 4-site 8-bath systems for the tt-Hux 
Hubbard model. In the honeycomb lattice, we find that 
even if we add the bath sites, the gap Agp opens at any 
infinitesimal U value when we use the 6-site hexagonal 
ring cluster as the reference system.In the TT-flux 
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Fig. 6. (Color online) Calculated single-particle gap Asp in the 
PM state of the (a) honeycomb and (b) tt-Aux Hubbard models. 
For the honeycomb lattice, we use the hexagonal 6-site cluster (h6) 
and 6-site 6-bath cluster (h6-6b). For the tt-Aux lattice, we use the 
square 4-site cluster (7r4), 4-site 4-bath cluster (7r4-4b), and 4-site 
8-bath cluster (7r4-8b). (c) tt-Aux lattice with renormalized hopping 
parameter violating the original translational symmetry and (d) 
its noninteracting single-particle gap as a function of 


Hubbard model, we also find that the gap Agp opens at 
any infinitesimal U value when we use the 4-site square 
cluster as the reference system. Therefore, even though 
we use two bath sites per correlation site, the gap opens 
at infinitesimal U values, which does not agree with the 
argument in Ref. 21 that at least two bath sites per cor¬ 
relation site are necessary to discuss the MIT in the hon¬ 
eycomb lattice. 

Rather, our results agree with the statement in Ref. 22 
that the opening of the gap at infinitesimal U/t values 
is not caused by the bath degrees of freedom but by 
the cluster geometry, which violates the original transla¬ 
tional symmetry of the lattice. We show the latter case 
in Fig. 6(c), where the original translational symmetry of 
the TT-flux lattice is violated by the renormalization of the 
hopping parameter t* by the interaction only within the 
cluster, leading to t* ^ t. Then, as shown in Fig. 6(d), 
the noninteracting band with t* and t has a finite single¬ 
particle gap unless =t. A similar discussion has been 
given for the honeycomb lattice,where the 6-site 
hexagonal clusters with the renormalized hopping pa¬ 
rameter are connected with the bare hopping param¬ 
eter t. A “plaquette insulator” state is thus realized at 
t* 7^ t in the noninteracting limit. This is the reason why 
the single-particle gap opens at infinitesimal U/t values. 
However, we here point out that it is always possible to 
make an appropriate choice of clusters that maintains 


the Dirac zero-gap situation even though it violates the 
original translational symmetry, examples of which are 
shown in Figs. 2(d)-2(f) where the gap does not open at 
small values of U. Thus, the statement in Ref. 22 is too 
strict. Careful choice of the clusters in the quantum clus¬ 
ter methods such as CDMFT, the VGA, and the CDIA 
enables one to discuss the MIT of Dirac fermion systems 
without spurious opening of the gap. 

4. Summary 

We have made a comparative study on the MIT of 
Dirac electrons in the honeycomb and tt-Aux Hubbard 
models using the VGA and CDIA, where we have calcu¬ 
lated the single-particle gap and staggered magnetization 
as functions of the interaction strength U. We have paid 
particular attention to the choice of the cluster geometry 
and the inclusion of the bath sites. We have thus con¬ 
firmed that the spurious single-particle gap that opens 
at infinitesimal U values is not caused by the bath de¬ 
grees of freedom but rather by the cluster geometry. We 
have shown that with increasing I/, the first-order MIT 
to the nonmagnetic MI phase occurs in the presence of 
electron-number fluctuation. However, the AFMI phase 
always preempts this MIT, at least in the present mod¬ 
els, and therefore the SL phase previously suggested to 
emerge between the Dirac SM and AFMI phases is ab¬ 
sent in these models. 

Our results imply that, if the AF ordering can be sup¬ 
pressed by, for example, to the effect of spin frustrations 
in the triangular and related lattices, one may expect 
that the MI phase without AF orders will preempt the 
AFMI phase, resulting in the emergence of the SL phase 
in the intermediate coupling region, as was pointed out 
recently in Ref. 35 for the triangular 7r-ffux Hubbard 
model. 

We thank K. Seki for enlightening discussions. T. K. 
acknowledges support from a JSPS Research Fellowship 
for Young Scientists. This work was supported in part by 
KAKENHI Grant No. 26400349 from JSPS of Japan. 


1) N. F. Mott, Metal-Insulator Transitions^ (Taylor & Francis, 
London,1990) 2nd ed. 

2) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 
1039 (1998). 

3) Note that the honeycomb Hubbard model we use is not suAi- 
cient for a realistic description of graphene although it is the 
simplest for reproducing its Dirac-type electrons: see, for ex¬ 
ample, A. Castro Neto, F. Guinea, N. Peres, K. Novoselov, 
and A. Geim, Rev. Mod. Phys. 81, 109 (2009); V. N. Ko¬ 
tov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro 
Neto, Rev. Mod. Phys. 84, 1067 (2012); M. Schiller, M. Rdsner, 
T. O. Wehling, A. 1. Lichtenstein, and M. 1. Katsnelson, Phys. 
Rev. Lett. Ill, 036601 (2013). 

4) S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992). 

5) 1. F. Herbut, Phys. Rev. Lett. 97, 146401 (2006). 

6) M.-T. Tran and K. Kuroki, Phys. Rev. B 79, 125125 (2009). 


5 





















J. Phys. Soc. Jpn. 


7) S. A. Jafari, Eur. Phys. J. B 68, 537 (2009). 

8) Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Mu- 
ramatsu. Nature 464, 847 (2010). 

9) S. Sorella, Y. Otsuka, and S. Yunoki, Sci. Rep. 2 , 992 (2012). 

10) F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013). 

11) F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. 
Herbut, arXiv:1411.2502. 

12) C.-C. Chang and R. Scalettar, Phys. Rev. Lett. 109, 026404 
(2012). 

13) D. Ixert, F. F. Assaad, and K. P. Schmidt, Phys. Rev. B 90, 
195133 (2014). 

14) W. Wu, Y.-H. Chen, H.-S. Tao, N.-H. Tong, and W.-M. Liu, 
Phys. Rev. B 82, 245102 (2010). 

15) A. Liebsch, Phys. Rev. B 83, 035113 (2011). 

16) A. Liebsch and H. Ishida, J. Phys.: Condens. Matter 24, 
053201 (2012). 

17) S.-L. Yu, X. Xie, and J.-X. Li, Phys. Rev. Lett. 107, 010401 

( 2011 ). 

18) W. Wu, S. Rachel, W.-M. Liu, and K. Le Hur, Phys. Rev. B 
85, 205102 (2012). 

19) R.-Q. He and Z.-Y. Lu, Phys. Rev. B 86, 045105 (2012). 

20) K. Seki and Y. Ohta, arXiv:1209.2101. 

21) S. Hassan and D. Senechal, Phys. Rev. Lett. 110, 096402 
(2013). 

22) A. Liebsch and W. Wu, Phys. Rev. B 87, 205127 (2013). 


23) Q. Chen, C. H. Booth, S. Sharma, C. Knizia, and C. K.-L. 
Chan, Phys. Rev. B 89, 165134 (2014). 

24) M. Laubach, J. Reuther, R. Thomale, and S. Rachel, Phys. 
Rev. B 90, 165136 (2014). 

25) Y. Hatsugai, T. Fukui, and H. Aoki, Phys. Rev. B 74, 205414 
(2006). 

26) M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. 
91, 206402 (2003). 

27) D. Senechal, arXiv:0806.2690. 

28) M. Potthoff, in Strongly Correlated Systems - Theoretieal 
Methods, ed. A. Avella and F. Mancini (Springer, Berlin, 2012) 
Vol. 171, Chap. 10, p. 303. 

29) M. Potthoff, Eur. Phys. J. B 32, 429 (2003). 

30) M. Potthoff, Eur. Phys. J. B 36, 335 (2003). 

31) C. Dahnken, M. Aichhorn, W. Hanke, E. Arrigoni, and 
M. Potthoff, Phys. Rev. B 70, 245110 (2004). 

32) D. Senechal, in Strongly Correlated Systems - Theoretieal 
Methods, ed. A. Avella and F. Mancini (Springer, Berlin, 2012) 
Vol. 171, Chap. 11, p. 341. 

33) M. Balzer, B. Kyung, D. Senechal, A.-M. S. Tremblay, and 
M. Potthoff, Europhys. Lett. 85, 17002 (2009). 

34) D. Senechal, D. Perez, and M. Pioro-Ladriere, Phys. Rev. Lett. 
84, 522 (2000). 

35) S. Rachel, M. Laubach, J. Reuther, and R. Thomale, 
arXiv:1411.4649. 



